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A passive scalar is advected by a velocity field, with a nonuniform spatial source that 
maintains concentration inhomogeneities. For example, the scalar could be temperature with 
a source consisting of hot and cold spots, such that the mean temperature is constant. Which 
source distributions are best mixed by this velocity field? This question has a straightforward 
yet rich answer that is relevant to real mixing problems. We use a multiscale measure of 
steady-state enhancement to mixing and optimize it by a variational approach. We then solve 
the resulting Euler-Lagrange equation for a perturbed uniform flow and for simple cellular 
flows. The optimal source distributions have many broad features that are as expected: they 
avoid stagnation points, favor regions of fast flow, and their contours are aligned such that 
the flow blows hot spots onto cold and vice versa. However, the detailed structure varies 
widely with diffusivity and other problem parameters. Though these are model problems, 
the optimization procedure is simple enough to be adapted to more complex situations. 



I. INTRODUCTION 

Consider a passive scalar advected by a velocity field, in the presence of an inhomogeneous spa- 
tial source. An obvious question is, which velocity fields are best at homogenizing the concentration 
field? This is a challenging question, and here we turn it around into a less obvious one: given a 
velocity field, which source distributions are best mixed by this field? For example, in a room with 
a given airflow, where to best position heating units to achieve a homogeneous temperature? We 
will see that answering this type of question gives considerable insight into optimal stirring flows 
in general. 

To carry out the optimization we need to narrow the problem further. First, we fix the ampli- 
tude of the velocity field, such as by specifying its total energy (in a bounded domain) or energy 
density (in an unbounded domain). Second, we restrict the source and velocity fields to be time- 
independent, so that the problem is steady. Finally, and most importantly, we need to specify how 
we measure the mixing enhancement due to stirring. Here we use a generalization of the variance 
of the concentration field, which has long been a popular measure of mixing. The reasoning is that 
a velocity field that is efficient at stirring should suppress fluctuations in the concentration field, 
and the variance decreases as these fluctuations become smaller. 

Previous work on optimization of mixing has focused on breaking flow symmetries [1-4], or 
optimizing quantities associated with chaotic advection, such as Lyapunov exponents or topological 
entropy [5-10]. Recent work has also involved optimizing the norm of the concentration field [11, 
12], in a manner similar to here. However, all these approaches differ from our work in that they 
do not involve body sources of scalar continually replenishing the variance of the concentration 
field. The optimal solutions we find are quite different, especially in that they do not tend to lead 
to creation of small spatial scales in the concentration field, but rather magnify the importance of 
efficient transport of temperature from sources to sinks. 

Motivated by Doering and Thiffeault [13] and Shaw et al. [14], we introduce a one-parameter 
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family of measures of mixing. The mixing enhancement factor £ p is thus defined by 

,_ ||(-A)^|| 2 
P ' ll(-A)^|| 2 ' Uj 

where 0(x) is the concentration of the advected scalar, ||-|| 2 is the L 2 norm, and A the Laplacian. 
The 9{x) in the numerator is the reference concentration obtained for the same source distribution 
and diffusivity, but in the absence of stirring. The enhancement factor thus tells us how much 
better the velocity field is at suppressing fluctuations than if we didn't stir at all. The above 
definition is appropriate for the steady advection-diffusion problem that we consider in this paper. 
The relevant definition for the time-dependent problem is given in Doering and Thiffeault [13], 
and involves taking the long-time average of the numerator and denominator in (1). As mentioned 
above, for simplicity we shall restrict our study to time-independent flows and sources, though 
the more general formulation is not conceptually more complicated. We assume without loss of 
generality that the source and initial condition, and consequently the scalar concentration, have 
spatial-mean zero. For concreteness, we will often refer to 9 as 'temperature' or 'heat', and speak 
of 'hot' and 'cold' regions, but the considerations here apply to any passive scalar. 

The numerator in the definition of the enhancement (1) is designed to avoid pathological solu- 
tions, such as ones that concentrate the source at very small scales. Nevertheless, our definition (1) 
sometimes behaves counterintuitively: for instance, the enhancement can be increased by rearrang- 
ing the source on large scales to penalize diffusion. These cases are not the most relevant ones, and 
do not justify introducing a new definition of the enhancement from that used previously Doering 
and Thiffeault [13], Shaw et al. [14], Thiffeault et al. [15]. Moreover, the definition (1) is well-suited 
to optimizing the velocity field for a fixed source distribution, which is a more important (and much 
more difficult) problem for which partial results are available. Keeping the same definition allows 
comparison to these earlier results. 

For p = 0, the enhancement factor (1) involves the scalar variance and was introduced by Thif- 
feault et al. [15]. Varying p preferentially weighs the smaller (p > 0) or larger (p < 0) scales, 
providing a different measure of mixing enhancement. This measure of mixing efficiency was used 
by Doering and Thiffeault [13] and Shaw et al. [14] in the context of statistically-steady turbulent 
flows, and they found that the scaling of the enhancement factor with the energy of the velocity 
field depends strongly on the nature of the source. For p = — 1, the enhancement factor is closely 
related to the mix-norm [11, 12]. 

In this paper we maximize (1) using a variational approach (Section II). The variation leads to 
an Euler-Lagrange equation where the leading eigenvalue is the optimal enhancement factor, with 
the corresponding eigenfunction giving the optimal source distribution. Some salient features of 
the optimal solutions are that (i) the optimal source distribution avoids having hot or cold spots 
over stagnation points of the flow; (ii) regions of high velocity are favored; and (iii) the hot and cold 
spots are positioned such that the flow sweeps hot onto cold and vice versa. Though most solutions 
have these broad features, they differ considerably in their details and often change dramatically 
(but continuously) as parameters, such as the diffusivity or the exponent p in (1), are varied. 

We illustrate the range of solutions by considering first the simplest situation, that of a uniform 
velocity field, as discussed in Plasting and Young [16], Doering and Thiffeault [13], and Shaw et al. 
[14]. In that case the optimization problem can be solved analytically (Section III). We then 
impose a perturbation to the uniform flow, leading to either a shear flow or a wavy flow, and solve 
for the optimal source using perturbation theory. The resulting optimal sources favor regions of 
high velocity in the shear flow, but the wavy flow leaves the optimal solution unchanged from the 
uniform flow, and only decreases its efficiency. This shows, unsurprisingly, that it is still possible 
to improve the enhancement factor after optimizing the source by optimizing the velocity field. 
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We then move to direct numerical solution of the optimization problem for a simple cellular 
flow (Section IV). The basic problem turns out to be doubly-degenerate in that there are two 
independent source distributions that give the same optimal enhancement. This is a consequence of 
a symmetry of the flow, and we verify that breaking this symmetry by adding a small perturbation 
to the velocity field selects a unique optimal solution. We also use the cellular flow to study 
the range of behavior as the diffusivity and the exponent p in (1) are varied. In all these cases, 
the optimal source distribution converges at extreme parameter values (large or small diffusivity 
or p — > ±00). For all the cases, we also compare the optimized enhancement factor to two simple 
reference sources, sinx and cosx, chosen to capture the spatial-phase dependence of the optimal 
source distribution. Because of the phase of the rolls, the cosx is much more efficient than sinx, 
and its enhancement factor is very close to optimal. This is not in itself a drawback, since it shows 
that our definition of enhancement is fairly robust to changes in the source, a desirable feature for 
practical implementation. Finally, we offer some concluding remarks in Section V. 



II. THE OPTIMIZATION PROCEDURE 



We consider the time-independent advection-diffusion equation for a passive scalar with con- 
centration 9(x), spatial source s(x), and diffusivity k, 

u(x) -V0 - kA6 = s(x), (2) 

in [0, L] d with periodic boundary conditions. The (given) velocity field u(x) is incompressible, 
V • u = 0. Both u{x) and s{x) are assumed to be sufficiently smooth. We assume that the source 
and initial condition have spatial mean zero, which implies that the scalar concentration also has 
mean zero. We remark that the solution T(x,t) of the evolution problem 

d t T(x, t) + u(x) ■ VT(x, t) - kAT(x, t) = s{x), 

converges, in the limit t — ► +00, to the solution 9(x) of the steady problem (2), the convergence 
being strong in L 2 . Hence, for steady sources and stirrers, it is sufficient to consider the stationary 
problem (2). 



A. Optimal Mixing Enhancement 

Our goal is to maximize the enhancement factor £ p defined by (1), 

P ' IK-a^H 2 .' 

where 9 solves equation (2) in the absence of advection, 

-kA6 = s, (3) 

with periodic boundary conditions on [0,L] rf . We assume that the velocity field is given. In 
maximizing the enhancement factor £ p , we fix the L 2 norm of the source and of the velocity field 
(or equivalently, the kinetic energy of the flow). 
Define the linear operators 

L := u(x) ■ V — kA and £ := — kA, 
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which implies 



from which we can write the solution to (2) and (3) as 

6 = £- 1 s and 9 = Z- 1 s. 
We can then rewrite the enhancement factor (1) as 

P \\(-A)p/^s\\ 2 2 (sA^s)' U 

where the self-adjoint operators A p and Ap are 

Ap := JL(-A)- p L* , A p :=Z{-A)- p Z* = k 2 {-A) 2 -' p , (5) 

and we have used the notation (•) to denote integration over [0, L] d . To maximize £ 2 , we compute 
its variation with respect to s and set it equal to zero, 

K = " £ p V s ) 6s ) = °> ( 6 ) 

Ap 1 s = E-pAp 1 s , (7) 

t/lp./lp 5 = £p 5. (^) 

This is an eigenvalue problem for the operator X p := A V A V X . The optimal source is given by 
the ground state of the inverse of this operator, and the normalized variance is given by the 
corresponding (first) eigenvalue. 

The operators Ap 1 and Ap 1 are self-adjoint from L 2 ([0, L] d ) to L 2 ([0, L] d ); furthermore, they 
are both positive operators in L 2 ([0,L] d ) (restricted to functions with mean zero). Consequently, 
the generalized eigenvalue problem (7) has real positive eigenvalues, and the eigenfunctions s 
and s' corresponding to distinct eigenvalues are orthogonal with respect to the weighted inner 
product (s , s') := (sA^s'}. 

Our goal now is to calculate the optimal source and the corresponding mixing enhancement 
factor for some simple velocity fields. In particular, in Section III we will consider a weakly 
perturbed uniform flow in two dimensions, and in Section IV we will consider cellular flows. We 
will be primarily concerned with the eigenvalue problem (7) or (8) for p = 0, i.e. 

AA.- 1 s = £. 2 s, (9) 

with A := .Ao = A := Ao = n 2 A 2 and £ := £o- Notice that the operator A^ 1 is a diagonal 
operator in Fourier space with entries n 2 \k\~ 4 , where k is the wavevector. The large negative power 
of \k\ indicates that this operator acts as a low-pass filter, suppressing high frequencies. We shall 
return to the case p / in Section IV C. 



or 



B. Is the Enhancement Maximal? 



Before considering specific examples in Sections III and IV, let us demonstrate that the optimal 
solution obtained in Section II A is indeed a global maximum. We could do this using the second 
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variation of £p, but instead we proceed more directly by expanding an arbitrary source in terms 
of eigenfunctions satisfying Eq. (7), 

s(x) = J2a (i) s (i \x), (10) 

i 

and inserting into expression (4) for the mixing enhancement factor, 



(11) 



where we used the eigenfunction property. Now use the orthogonality property of the 
(sW Ap 1 s(fi) = Sij, and the fact that £ p > tp\ that is the optimal enhancement factor £ p is 
the largest eigenvalue, to find 

<fj^f) = Ej( Q(i) ) 2 ( g(i) ^p lg(i) ) Ej(« (i) ) 2 ( g(i) -^ P lg(i) ) _ £2 

which proves that £ p is indeed optimal, since the enhancement factor of no other source can exceed 
it. This is a global argument, which relies only on the eigenvalue problem (7) providing a complete 
set of orthogonal eigenfunctions. 



III. UNIFORM FLOW WITH PERTURBATION 



As a first test case for the optimization procedure of Section II, we consider a uniform flow (con- 
stant magnitude and direction in space). The uniform flow illustrates a fundamental mechanism 
involved in source optimization, as mentioned in Plasting and Young [16] and Shaw et al. [14]: the 
velocity field sweeps the hot source onto the cold sink, and vice versa. More generally, we expect 
that the optimal source will tend to have contours perpendicular to the flow. The uniform flow 
maximizes the mixing enhancement factor for a one-dimensional source, at fixed kinetic energy. 

To capture another important feature of optimal sources, in Section III B we perturb the uniform 
flow to make either a shear flow or a wavy flow. The shear flow perturbation will show that optimal 
sources are localized over rapid regions of the flow. The wavy flow perturbation will show that 
sometimes aligning the source contours perpendicular to the flow is too 'costly', and the optimal 
solution is left unchanged from the uniform flow. The cost incurred is due to the weighting of the 
enhancement factor (1) by the purely-diffusive solution: a change in the source might make it more 
efficient, but it can also make the purely-diffusive solution more efficient, yielding no net gain. 

Note that in this section we will restrict our optimization to the variance, that is with p = in 
the mixing enhancement (1). We shall return to the effect of varying p in Section IV C. 



A. Uniform Flow 



For completeness, in this short section we essentially rederive the results of Plasting and Young 
[16], Doering and Thiffeault [13], and Shaw et al. [14] on the optimality of a uniform flow. We 
start with a time- independent spatially uniform flow on [0,L] d , 



u(x) = U. 
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From definition (5) with p = 0, we have that 

A = k 2 A 2 — (U ■ V) 2 , A = k 2 A 2 , 

and 

X = AAT 1 = / - k- 2 (U ■ V) 2 A~ 2 , 

Since DC is a differential operator with constant coefficients on a periodic domain, its eigenfunctions 
are easily verified to be 

Sfc(aj) = W 8 ' 35 + c.c. , 

where k is the wavevector and s& is a normalization constant, and the corresponding eigenvalues 
are 

\U-k\ 2 

+ K 2 |fc| 4 • 

We can maximize the above expression by choosing k to be the shortest allowable vector parallel 
to the uniform flow U. (Not all magnitudes and directions are allowed for wavevectors because the 
domain is periodic: the components of k are integer multiples of 2n/L.) 

In the particular case where the uniform flow U is along the x-axis, U = U e x , we have that 

I TJ 2 L 2 

where we have defined the Peclet number Pe := UL/2ttk. As expected, in the limit as k — ► oo 
(Pe — > 0), the mixing enhancement factor converges to 1, i.e. in the large diffusivity limit the 
enhancement due to a flow becomes negligible. On the other hand, in the limit n — > (Pe — > oo) 
the enhancement factor grows like 



B. Shear and Wavy Flows 

Consider now a two-dimensional uniform flow along the x-axis perturbed by a weak flow, 

|| Wl ||2 

u(x,y) = Ue x + eu 1 (x,y) - e 2 ^j^e x , (13) 

with e<Cl and (tti) = 0. The second-order term is included to make ||tt|| 2 = U + 0(e 4 ), so that 
we can compare the effect of the perturbed and unperturbed flows at equal amplitude. Specifically, 
we consider perturbations of the form 

m(x, y) = u lx (y) e x + u ly (x) e y , 

which are simple to analyze but yield important insight. Because the base flow is in the e x direction, 
the u\ x {y)e x term is a shear flow perturbation, and the u\ y {x)e y term is a wavy flow perturbation. 
Our goal is to obtain an asymptotic expansion for the optimal source and enhancement factor. For 
this we need to study perturbatively the eigenvalue problem (9), which we write as 



J\.J\ s — A s ; 



(14) 
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with A = A = £X* = k 2 A 2 , and A = £ 2 . We have that 

& = £jq + e£j\ + e 2 /L2 , 

with 

£<o = Ud x - kA , Li = ui ■ V , £ 2 = u 2x d x , 
where u 2x is the coefficient of the second-order term in (13). Consequently, 

A = ^o^o + ^(^i^o ~~ ^o^i) + £ 2 (^2^o ~~ ^0^2 — ^l^l) 
=: Aq + eAi + e 2 A 2 . 

Note that in this section the subscripts on the A's correspond to their order in e, and not to the p 
subscript as in Eq. (5) (recall that p = for the present section). We expand s and A in a power 
series in e, 

s = s + es\ + e 2 s 2 + . . . , 
A = A + e\i + e 2 X 2 + 

insert the expansions for A, s and A in equation (14), and equate like powers of e to obtain the 
sequence of equations 

A A~ 1 s = \ s , (15a) 
AtA^so + A A~ 1 s 1 = A si + Ais , (15b) 
A 2 A~ 1 sq + AiA _1 si + AqA~ 1 s 2 = Ao«2 + Ai«i + \ 2 sq . (15c) 

Consider first equation (15a), the unperturbed equation: from Section III A we have that 

U 2 

so = s e lk ° x + c.c. , A = 1 + = 1 + Pe 2 . 

where k$ := 2n/L. We proceed now with equation (15b): multiply by s$ and integrate over [0, L] 2 
to obtain 

(so&xHqA^so) - (so Lo&iA^ 1 s ) + (soAoA^si) = A (s si) + Ax || s 111 • (16) 

Notice that Ao is a differential operator with constant coefficients, and hence it commutes with 
A^ 1 . Furthermore, both operators are self-adjoint, so that 

(s AqA^sx) = (sqA^Aosi) = ((AqA^so) si) = A (s si) (17) 
and equation (16) simplifies to 

(so^i^o-^ -1 ^) - (s £ ^i-A _1 so) = AiHsolla 

leading to 

(so^i^o-^ -1 ^) {so&o&iA^so) 

1 ~ ii — ii2 ii — ii2 ~~ ' 

II s o II2 \\ s o\\ 2 



8 



where the last equality follows from a straightforward calculation. This is a consequence of the 
invariance of the enhancement factor under reversal of the pertubation, from which all odd powers 
of e in the enhancement factor will vanish. 

Now we compute s\. We set Ai = in equation (15b) to obtain 

(A - AqA' 1 ) si = AiA _1 s - (18) 
For definiteness, we specify the form of the perturbation, 

ui x (y) = in x e ikiy + c.c. 

Notice that we do not yet need to specify u\ y (x): it only affects the result at the next order. We 
look for a solution to (18) in the form of 



8l (x,y) = A k0 k,hui x e i{k0X+kiy) + B kokl s u lx e^ x ~ k ^ + c.c, (19) 



and find 



Let us proceed now with the calculation of A2. We multiply (15c) by so, integrate over [0, L] 2 
and use (17) to obtain 

. (s A 2 A- 1 so) , (so^-^i) 

A2 ~~ m [72 1 n jT2 ' 

ll s o|| 2 INH2 

which after substituting our solution (19)-(20) for s\ simplifies to 

(Ak 2 U 2 + kjK 2 ) , |2 1 „ ||2 

A2 " U 2 K 2 k 2 (2k 2 + k 2 ) llUlxh n 2 k 2 llUlyh ' 



We are really after the enhancement factor, which is 



£ = VAo + e 2 A 2 + 0(e 4 ) = + e 2 -^= + 0(e 4 ) , 



so after putting it all together, we can write £ in terms of dimensionless quantities, 

F-,/TT^ + S / + « 4 Pe- 2 ) \Mj _ JKJ||1 Q(g4) (21) 

where a '■= k±/ko and recall that Pe = U/nko. 

Some observations regarding (21) are in order. First, notice that a nonzero u\ x (shear flow 
perturbation) increases the enhancement factor, while a nonzero u\ y (wavy flow) decreases it. Sec- 
ond, the form of u\ y is irrelevant, and it does not affect the first-order optimal source distribution. 
Third, in the limit of large k (small Pe), the Pe -2 term in (21), arising from k 2 A 2 in Aq, invalidates 
the expansion. We discuss each type of perturbation in turn. 

Figure 1(a) shows the optimal source distribution for U = 1, L = 2ir, k = 0.01 (Pe = 100), 
u lx = l/s/2, k\ = 2, e = 0.15, u\ y = 0. This is a shear flow, with contour lines shown in the 
background. The contour lines are closer together at the midpoint, indicating faster flow, since 
the perturbation is u\ x {y) = \/2cos2y. Accordingly, the optimal source is localized at these points 
of faster flow. However, there is a limit to how localized it can get: if the source were to bunch 
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(a) (b) 

FIG. 1: Optimal source distribution for the flow e x + e u±(x, y) — e 2 \ \\u\ || 2 e x with e = 0.15, k = 0.01 and 
(a) Ui(y) = \/2cos2y e x ; (b) U\{x) = \/2 cos 2x e y . The horizontal phase of the solution is arbitrary. The 
background shading shows hot (red, or dark gray) and cold (blue, or light gray) regions, separated by tepid 
regions (white). The contour lines are streamlines of the flow. 

up too much in the faster region, the purely-diffusive solution would be more effective (i.e., have 
lower variance), reducing the enhancement factor. As the perturbation gets larger, the optimal 
source will bunch up more in the faster regions, since the gain becomes greater. The enhancement 
factor is always greater than in the unperturbed case, suggesting that the advantage of faster flow 
regions outweighs the disadvantage of localizing the source. Finally, note that for large a (large 
perturbation wavenumber k\), the mixing enhancement factor becomes independent of a, reflecting 
the fact that a fine-scale perturbation is overwhelmed by diffusion. 

Conversely, Figure 1(b) shows the optimal source distribution for the same parameters as the 
shear flow above but with perturbation u\ y = \/2cos2x, u\ x = 0. The streamlines in the back- 
ground indicate that this is a wavy flow, with no dependence on y. The optimal source distribution, 
again shown in the background (here as cosj;, but the phase is arbitrary), is exactly the same as 
in the absence of perturbation. Intuitively, one could expect the source to be more efficient if it 
tilted to present its contours perpendicular to the flow, giving a wavy source. This is not the case, 
since the purely-diffusive solution would be more effective for a wavy source, and thus lower the 
enhancement factor since it enters the numerator in the definition (1). Localizing the source in 
regions of faster flow (as for the perturbation in Fig. 1) is a more important effect than aligning 
the contours of the source perpendicular to the flow. 

IV. NUMERICAL RESULTS FOR THE CELLULAR FLOW 

In Section III we derived expressions for optimal sources and mixing enhancement factors for 
perturbations of a uniform flow. More general velocity fields require a numerical approach. To 
maximize the enhancement factor, we again have to solve the generalized eigenvalue problem (8) 
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for the optimal source distribution s. For numerical implementation, it is preferable to solve the 
equivalent self-adjoint eigenvalue problem 

(A- 1 ' 2 AA- 1 ' 2 )r = £ 2 r, s = A 1 ' 2 r, (22) 

for the eigenvector r, which then yields the optimal source distribution s. The advantage of the 
form (22) is that the self-adjoint structure of the operator is explicit. In practice, we expand r as a 
Fourier series, so that A is a diagonal matrix. We then solve (22) using Matlab's eigs routine for 
sparse matrices. We set the box size L = 2ir throughout this section. Note that we shall deal only 
with p = (variance optimization) until Section IV C, so we leave off the p subscript until then. 

A. Cellular Flow 

We consider the perturbed cellular flow with streamfunction 

V>(x, y) = A (sin x sin y + 5\ sin 2x + 82 sin 2x sin 2y) (23) 

with velocity field u = (u x ,u y ) = (d y ip,—d x ip), where A is a normalization constant chosen to 
make ||w||2 = 1- With 8± = 82 = (the basic cellular flow), the operator £ is invariant under the 
discrete symmetry group S generated by the transformations 

G\ (x,y) = {y , —x + ir), rotation with vertical translation; (24a) 

G2 (a;, y) = (x + 7r , y + n), diagonal translation. (24b) 

The Abelian group S has order 8 and is characterized by G\ = G\ = I, G1G2 = G2G1. 1 For 
either 8\ or 82 nonzero, the perturbations break the symmetry G\ of the cellular flow, and we shall 
look at their effect in turn. 

First we set 8\ = 82 = in (23) and solve (22) with k = 0.01. In this case, there are two 
independent optimal source eigenfunctions with degenerate enhancement factor £ = 87.34, shown 
in Fig. 2. In the foreground are contour lines of the streamfunction. Any superposition of the two 
eigenfunctions in Fig. 2 will give the same enhancement factor. The degeneracy is a consequence 
of the symmetry G\ of the cellular flow: indeed, the two eigenfunctions are related to each other 
(up to a sign) by the transformation (24a). The two eigenfunctions also separately have the G2 
symmetry. This situation, where the eigenfunctions corresponding to the same eigenvalue are 
related by a unitary representation of the symmetry group of an operator, is familiar from quantum 
mechanics [17, 18]. 

For comparison, the enhancement factor for the same flow but with the reference source s(x) = 
since is 50.01, and with the reference source cosx is 86.61. Hence, for sinx the optimal source gives 
a 74.7% improvement in the enhancement factor, but only 0.9% for cosx. It is remarkable how 
close cos x comes to the optimal enhancement factor, which shows that the optimal source is rather 
'robust', so that small changes in its shape do not lead to huge changes in the enhancement factor. 
From a design standpoint, this is highly desirable. Table I summarizes the optimal enhancement 
factor results for different values of the perturbations 8\ and 82- 

How to interpret the optimal source distribution in Fig. 2? The lesson learned from the uniform 
flow of Section III is that, ideally, the source should be such that the velocity field advects heat 
from hot to cold and vice versa. This is clearly happening in Fig. 2 to some extent, since in both 
eigenfunctions the hot and cold spots are situated such that the rolls easily advect heat between 



1 S is the direct product of a cyclic group of order 4 and a cyclic group of order 2. 
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s (l) s (2) 




FIG. 2: Optimal source distributions for the pure cellular flow (23) with Si = 62 = 0. The two pictures show 
degenerate orthogonal eigcnfunctions with enhancement factor £ = 87.34. Note how there is no source of 
heat over the stagnation points. (See the caption to Fig. 1 for a key to the background shading.) 

TABLE I: Optimal mixing enhancement factors £ ptimai with p = for the perturbed cellular flow (23), 
with k = 0.01. Compare to the reference enhancements for a sinx and cos a; source distribution: the number 
in parentheses is the % improvement of the optimal source. The cosx source always does much better 
than sin a; because it straddles the rolls, whereas sin a; has hot and cold segregated into different rolls. 



Si 


§2 


^optimal 


£ si „ (%) 


£cos (%) 


note 








87.34 


50.01 (74.7%) 


86.61 (0.9%) 


Fig. 2 (degenerate) 


0.05 





87.59 


49.76 (76.0%) 


86.18 (1.6%) 


Fig. 4(a) 





0.05 


90.10 


50.26 (79.3%) 


86.47 (4.2%) 


Fig. 4(b) 


0.05 


0.05 


90.19 


50.01 (80.3%) 


86.04 (4.8%) 




0.2 





90.43 


46.43 (94.7%) 


80.41 (12.5%) 







0.2 


95.34 


53.35 (78.7%) 


84.59 (12.7%) 




0.2 


0.2 


97.60 


50.01 (95.2%) 


79.30 (23.1%) 


Fig. 4(c) 


0.5 





95.99 


35.37 (171%) 


61.25 (56.7%) 


Fig. 4(d) 





0.5 


94.91 


61.25 (55.0%) 


79.06 (20.1%) 




0.5 


0.5 


106.8 


50.01 (114%) 


64.55 (65.4%) 





hot and cold. However, this is not the whole story: the optimal source is also distributed so as 
to take advantage of regions of faster flow. This is readily apparent in Fig. 3, which shows the 
magnitude of the velocity field in the background, and contours of the optimal source distribution 
from Fig. 2 (right). The hot and cold spots (elliptic regions) are clearly localized over the regions 
of rapid flow (pale background). Some of the fast regions appear to have no hot or cold spots, but 
these regions are favored in the other eigenfunction, so the symmetry is respected. 

Next, we resolve the degeneracy of the optimal source distribution by setting 5\ = 0.05 in (23) 
while keeping 62 = 0, giving the streamfunction shown as contour lines in the foreground of 
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FIG. 3: The background shows the magnitude \u\ of the velocity field for the pure cellular flow, Eq. (23) 
with 5± = 82 = 0. The contour lines are from the optimal source distribution on the right in Fig. 2. Note 
how the source-sink pairs (rolls in the contours) are clustered over regions of high speed (pale) and avoid 
the stagnation points (dark) . 

Fig. 4(a). The rolls now have a slight asymmetry that breaks the rotational symmetry G\ of 
the unperturbed cellular flow. As a consequence, the normalized optimal eigenfunction is now 
unique, and is shown as the shaded background in Fig. 4(a). It is very close to the degenerate 
eigenfunction on the right in Fig. 2, and converges to it as #1 — * 0. The enhancement factor in this 
case is almost unchanged, 87.59. Again, the cosine reference source has much higher enhancement 
factor (86.18, 1.6% improvement for the optimal source) than the sine reference source (49.76, 
76.0% improvement for the optimal source). 

Finally, we set 81 = and 62 = 0.05 in (23), with the streamfunction shown as contours in the 
foreground of Fig. 4(b). As for the previous perturbation, this one also breaks the G\ symmetry 
and causes the optimal eigenfunction to become unique, but this time a superposition of the two 
degenerate eigenfunctions in Fig. 2 is selected. The optimal enhancement factor, 90.10, is again 
almost unchanged by this small perturbation. 

To summarize this section, we presented three cases at fixed diffusivity k = 0.01. The first was 
the cellular flow, for which we get doubly-degenerate optimal eigenfunctions. We then presented 
two symmetry-breaking perturbations in turn, showing how these select a particular mixture of 
the degenerate eigenfunctions to create a unique optimal source distribution. Since the perturba- 
tions are small, the optimal enhancement factor in all these cases is about the same, showing an 
improvement of about 75-80% over the reference source sinx, but only 1-5% over the reference 
source cosx. This latter modest improvement is best seen not as a failure of the optimization 
procedure, but as an advantage, since robustness is always desirable. In fact robustness can easily 
be gauged by looking at the magnitude of the next largest eigenvalues, to see how far they are 
from the dominant one(s). Table I also shows that the modest improvement over the cosine source 
is an accident, since for many other velocity fields improvements well above 50% are seen for both 
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(c) (d) 

FIG. 4: Optimal source distribution for the perturbed cellular flow (23) with (a) 8i — 0.05, S 2 — 0, mixing 
enhancement factor £ = 87.59; (b) Si = 0, 5 2 = 0.05, £ = 90.10; (c) 8 1 = 5 2 = 0.2, £ = 97.60; (d) £1 = 0.5, 
<$2 = 0, £ = 95.99. The perturbations all break the G\ symmetry and selects a linear combination of the 
degenerate eigenfunction in Fig. 2. (See the caption to Fig. 1 for a key to the background shading.) 

sine and cosine. Two further examples for large perturbations are shown in Figs. 4(d) and 4(c). 

B. Dependence on Diffusivity 

Now we will fix 8\ = 82 = in (23) and vary the diffusivity, k. For k = 0.01, Fig. 2 shows 
the two degenerate optimal source eigenfunctions, and we will follow the change in the one on the 
left as k is varied. In Fig. 5 we show the change in the optimal source as k is increased from 0.1 
to 100. The optimal source distribution appears to become independent of k both for small k and 
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FIG. 5: For the same flow as in Fig. 2, optimal source distribution for different values of the diffusivity k, 
increasing from top left to bottom right. The cigcnfunction is doubly-dcgcncratc and corresponds to the one 
on the left in Fig. 2. For both small and large k the optimal source converges to an invariant eigenfunction. 
In all cases there is no source of temperature over the elliptic stagnation points, but in the large k case there 
are sources and sinks over some hyperbolic points. (See the caption to Fig. 1 for a key to the background 
shading.) 

large k, but the distributions are different. The transition between the two regimes occurs when k 
is of order unity. Though the two asymptotic sources are very different, they respect the general 
principles laid out in Section IV A: the source is arranged for effective transport of hot onto cold 
and vice versa, and regions of high speed are favored. In particular, note that the center of the 
rolls has a nearly zero, flat source distribution in all cases. 

Another perhaps surprising aspect of the large k solution in Fig. 5 is that it has complicated 
structure. In this large diffusivity limit, one would expect diffusion to dominate and gradients 
to be smoothed out. But since our mixing enhancement factor (1) compares the variance to the 
unstirred case, which already has very low variance, any amount of improvement will count. Hence, 
the complicated source for large k in Fig. 5 only gives a minute improvement to the enhancement 
factor. The large k optimal solution is particular in that it has some hot and cold spots localized 
over hyperbolic stagnation points. This is probably due to the high speeds along the separatrices 
being favored, even at the cost of straddling hyperbolic stagnation points a little. 

Figure 6 shows how the mixing enhancement factor varies as a function of the diffusivity. The 
solid line is for the optimal source, the dashed lines for the reference sources sin a; and cos x. For 
small k, the enhancement factor of all sources scales as this is the 'classical' scaling discussed 
in [13-15], where the enhancement factor is linear in the Peclet number. It has been rigorously 
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FIG. 6: For the flow with streamfunction as in Fig. 2, mixing enhancement factor £ — 1 as a function of the 
diffusivity k: optimal source (solid line), and sin a; and cos a; reference sources (dashed lines). For small k, 
the enhancement factor scales like For large k, the optimal solution returns to a k _1 approach to unity 
after a brief dip, while the reference source solution approaches unity as k~ 2 . 

proved in [13-15] that this scaling is optimal over all possible sources and velocity fields. 

For k near unity, the optimal enhancement factor has a small dip before converging towards 
unity as for large n. In contrast, the reference enhancement factors converges to unity as k~ 2 . 
This last scaling holds for the uniform flow of Section III when expanded in large k, and is verified 
for other flows and sources as well [13-15]. 

In summary, the optimal source distribution becomes independent of k for both large and small 
K, but of course for large k the efficiency gain is minimal (since the L 2 -norm of the velocity is fixed). 
For small k the efficiency gain is a constant multiple of the reference sources, but this multiple is 
fairly small for cosx (1.01), showing that optimization is very robust but not necessarily always 
worthwhile. Overall, the optimal enhancement factor scales as with a momentary break in 
the scaling that corresponds to the complicated change in topology seen in Fig. 5 for n near unity. 

C. Dependence on Exponent p 

Our final study will be to examine the behavior of the optimal enhancement factor as p is varied 
in (1). In Sections IVA-IVB we used p = 0; now we fix 5\ = 82 = 0, n = 0.01, and allow p to vary 
over negative and positive values. Figure 7 shows the optimal source distributions for p varying 
from —1 to 2. For both negative and positive p, the optimal source distribution converges rapidly to 
invariant patterns, and the two extremes in Fig. 7 are representative of those asymptotic patterns. 
The situation is thus entirely analogous to the case where diffusivity was varied (Fig. 5). 

The top-left picture in Fig. 7 (negative p) shows small, localized sources and sinks. In contrast, 
the bottom-right picture in Fig. 7 (positive p) shows large, regular localized sources and sinks. In 
fact, what is striking about the pattern is its simplicity: it is what one might take as a guess at 
an efficient source distribution, with no added frills. Thus, a high power of p might be useful in 
situations where a simple configuration is necessary due to engineering constraints. The reason 
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FIG. 7: For the same flow as in Fig. 2, optimal source distribution for different exponents p, increasing from 
top left to bottom right. The eigenfunction is doubly-degenerate and corresponds to the one on the left in 
Fig. 2. For both small and large p the optimal source converges to an invariant eigenfunction. In all cases 
there are no sources or sinks of temperature over the stagnation points. (See the caption to Fig. 1 for a key 
to the background shading.) 

for the simplicity is that spatial variations in the source favor the diffusion operator in £, and 
as p — * oo these are magnified. Thus, the source must remain as spatially simple as possible while 
trying to maximize alignment with the velocity. As p — > — oo, spatial variations of the source are 
downplayed by the norm, allowing more complexity. 

Figure 8 shows how the optimal mixing enhancement factor varies as a function of p. For p — > 
—oo, the enhancement factor goes to infinity, as does the enhancement factor of the two reference 
sources. For p — > oo, the enhancement factor also goes to infinity, but the reference source enhance- 
ment factors now approach constants (the constant is 1 for sinx). Again, we are seeing the effect 
of the diffusion term dominating when gradients are present, since these are amplified by (-A)p/ 2 . 
Note that the curve is symmetric about p = 1, which leads to a minimum there: whether this is 
true in general has not been determined, but we have not found a counterexample. In Appendix A 
we provide a partial proof by explicitly finding the symmetry between the operators Ai- V and A p , 
but only for large k. 
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FIG. 8: For the flow with streamfunction as in Fig. 2, mixing enhancement factor £ p - 1 as a function of 
the exponent p: optimal source (solid line), and sinx and cos a; reference sources (dashed line). The optimal 
enhancement factor is symmetric about p = 1, and for \p\ ^ 1 it grows as (2.2) > p >. 



V. DISCUSSION 

In both the perturbation problem (Section III) and the numerical examples (Section IV), the 
optimal source distributions tend to exhibit the following features: 

1. Avoidance of stagnation points of the flow, especially of elliptic type; 

2. Localization over regions of rapid flow; 

3. Alignment of the source contours perpendicular to the local velocity, so that hot is swept 
onto cold and vice versa. 

For a shear flow perturbation (Section IIIB), the optimal source bulges out over regions of faster 
flow. In contrast, a wavy flow perturbation leaves the optimal source unchanged over from a 
uniform flow, suggesting that localization is a more important effect than alignment. This also 
demonstrates that the optimal solution for source optimization does not necessarily correspond to 
the optimal solution for velocity optimization, given that optimal source. This is because at fixed 
energy the wavy flow always decreases the optimal enhancement factor from that of the uniform 
flow, for the same source distribution. 

These considerations show that the mixing enhancement factor and optimization procedure 
described in this paper (Section II) behave in a natural manner. Furthermore, the procedure yields 
a global maximum (Section II B), is numerically well-behaved, and is easy to implement. Indeed, a 
three-dimensional application to a real system is well within the realm of feasibility: the examples 
we presented here (except for the uniform flow) were two-dimensional, but for smooth flows the 
sparse nature of the advection-diffusion operator in Fourier space means that large problems can 
be solved with a modest amount of computer power and memory. 

The optimal source distribution becomes independent of the diffusivity k for both large and 
small k (Section IV B), and for exponent p [see Eq. (1)] negative and large or positive and large 
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(Section IV C). None of the source distributions achieved in these cases resemble each other. 
However, we observed that the eigenvalue spectrum of (7) is always symmetric about p = 1 (and 
has a minimum there) , which implies that the operators (and thus the optimal source distributions) 
are related by a unitary transformation we were unable to find in general (but see Appendix A 
for a partial result). The large positive p case is particularly interesting, since it favors source 
distributions that are very smooth, a desirable feature in practical applications. Another attractive 
feature we observed is the robustness of the optimal solution in the cases considered (Section IV A), 
but this will not necessarily hold in general. 

To widen the applicability of the procedure, a few complications will have to be introduced. 
First, time-dependence of the flow and the source is desirable, which will make the variational 
problem more difficult to solve in principle. Second, and more importantly, there remains the much 
more difficult problem of optimizing the velocity field given a source distribution. The variational 
problem is easily formulated, but does not present itself in the nice generalized eigenvalue problem 
[Eq. (7)] we saw here. The third and ultimate goal is a full dynamical coupling to the Navier-Stokes 
equation, with buoyancy and other effects included as appropriate. 

Throughout the paper we spoke of 'mixing' because our description involves the interplay of 
stirring and molecular diffusion. However, we highlighted the more interesting limit of small 
diffusivity k (large Peclet), since this is the situation in which stirring is more pertinent. In 
that case, it is clear that the pushing of hot fluid onto cold regions and vice versa is better 
described as 'transport' rather than 'mixing'. Thus, perhaps in that limit we should speak of a 
'transport enhancement factor', but since our analysis applies even for very large k we have kept 
the terminology, in spite of the fact that the velocity fields presented here do not lead to creation 
of small scales and thus may not 'mix' very well by some criteria. 
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APPENDIX A: SYMMETRY OF OPTIMAL MIXING ENHANCEMENT FACTOR 

ABOUT p = 1 

In this appendix we will motivate the symmetry of the optimal enhancement factor curve 
about p = 1, as evident in Fig. 8. We use the self-adjoint form (22): let 

%[u] :=A^ 2 A p [u]A^ 2 , 

where we have explicitly shown the dependence of the operators on u(x). We show below that 

9 B 2 -p[ti] IP = E p [3>u] + 0( K - 2 ) , (Al) 

where IP is the parity change operator, defined by r Pr{x) = r(—x) and "?u(x) = —u(—x), 
with T- 1 = IP. Since the spectrum is unchanged by the substitution u(x) — > — u(— x), estab- 
lishing (Al) proves that I^-pfw] and !B p [w] have the same spectrum, which explains the symmetry 
of the enhancement factor about p = 1. The eigenfunctions are related by a parity change. Unfor- 
tunately, Eq. (Al) is only an asymptotic result valid for large k, and we do not know the general 
form of the symmetry IP for smaller k, though numerical evidence suggests it exists. 
To show (Al) directly, we expand 

y^-pHx)} ■p = i+k~ 1 ((-A)-i ?u(x) ■ VIP (-A)f- 1 - (-A) 2- 1 y u (x) ■ vip (-A)-i) +0(k- 2 ) 

where we used the commutativity of A and IP. We also have 

%[-u{-x)} = 1 + k' 1 ((-A) _ 2 u(-aj) • V (-A)I- 1 - (-A)i" 1 u(-aj) • V (-A) _ i) + 0{k' 2 ). 

Now Eq. (Al) follows from !Pu(a;) • VIP = u(—x) ■ V, since both u and V reverse direction under 
parity change. 



